本章将对分析模拟轨迹的各种不同方法进行讨论, 并给出相应的分析程序的名称. 与这些程序相关的输入输出的具体信息, 可以利用www.gromacs.org的在线手册进行查阅. 输出文件通常都是Grace/Xmgr的图形格式.

首先, 在8.1节对分析中使用的组的概念进行解释. 8.1.2节将解释一个新的概念–动态选择–只有一些分析工具支持这一概念. 然后, 将对不同的分析工具进行讲解.

gmx make_ndx, gmx mk_angndx, gmx select

在第三章, 我们解释了如何在mdrun中使用 原子组 (参看3.3节). 对大多数分析程序, 你也必须选择原子组. 虽然大多数程序都能产生一些默认的索引组, 但这些程序总是可以从索引文件中读入组. 让我们考虑一个具体的例子, A和B双组分混合物的模拟. 当我们需要计算A对B的径向分布函数(RDF, Radial Distribution Function) \(g_{AB}(r)\) 时, 我们必须计算

\[4\pi r^2 g_{AB}(r) = V \sum_{i \in A}^{N_A} \sum_{j \in B}^{N_B} P(r) \tag{8.1}\]

其中, \(V\) 为体积, \(P(r)\) 为距离A原子 \(r\) 处发现B原子的概率.

通过让用户在一个简单的文件中定义组A和B的 原子编号, 我们就能以最通用的方式计算 \(g_{AB}\), 而无须在RDF程序中对粒子类型做任何假定.

因此, 组可以包含一系列 原子编号, 但在一些情况下也可以包含 分子编号. 我们还可以利用 三联原子编号 指定一系列键角, 利用 四联原子编号 指定一系列二面角, 利用 原子编号对 指定一系列键或(分子中的)矢量

需要对分析程序指定索引文件的类型. 许多程序都可帮你创建这样的索引文件(index.ndx), 它们在创建时使用你输入的构型或拓扑信息. 要创建包含一系列 原子编号 (如 \(g_{AB}\) 例子中)的索引文件, 可以使用 gmx make_ndx或gmx select. 要创建键角或二面角的索引文件, 可以使用gmx mk_angndx. 当然, 你也可以手动创建这些索引文件. 索引文件的通用格式如下:

可以看到, 组的名称写在方括号中. 随后的原子编号从1开始, 可以延续任意多行.

每个支持组的工具都会列出所有可选择的组, 供用户通过组的编号或名称进行选择. 实际上, 如果可以与其他的组区分开来, 使用组名称的头几个字母就足够了. 利用Unix shell可以在命令行中选择组名称, 而不是以通常的交互方式. 请查看www.gromacs.org给出的建议.

当没有为分析工具或grompp提供索引文件时, 会产生一些默认组以供选择:

System: 体系中所有的原子

Protein: 所有的蛋白质原子

Protein-H: 除氢原子外的所有蛋白质原子

C-alpha: \(\text C_\a\) 原子

Backbone: 蛋白质骨架原子, N, \(\text C_\a\) 和C

MainChain: 蛋白质主链上的原子: N, \(\text C_\a\), C和O, 包括C端的氧原子

MainChain+Cb: 蛋白质主链上的原子, 包含 \(\text C_\b\)

MainChain+H: 蛋白质主链上的原子, 包括骨架氨基氢原子和N端氢原子

SideChain: 蛋白质支链上的原子, 即除N, \(\text C_\a\), C, O, 骨架氨基氢原子, C端氧原子, N端氢原子以外的所有原子

SideChain-H: 除氢原子外的蛋白质支链上的原子

Prot_Masses: 除哑质量(创建NH3基团和色氨酸支链时使用的虚位点)外的蛋白质原子, 请参看5.2.2节. 只有当与Protein组不同时, 才包含此组

Non-Protein: 所有的非蛋白质原子

DNA: 所有的DNA原子

RNA: 所有的RNA原子

Water: 水分子(名称如SOL, WAT, HOH等). 参看residuetypes.dat文件中的完整列表

non-Water: 不属于Water组的任何粒子

Ion: 匹配residuetypes.dat文件中离子条目的任何名称

Water_and_Ions: Water和Ion组的组合

molecule_name: 所有不能被识别为蛋白质, DNA或RNA的残基/分子. 会为每个残基/分子创建一个组

Other: 不属于蛋白质, DNA, RNA的所有原子

不会创建空的组. 大部分组只包含蛋白质原子. 如果一个原子的残基名称存在于residuetypes.dat文件中, 并且被视为“蛋白质”, 这个原子就被视为蛋白质原子. 决定DNA原子, RNA原子等的方法类似. 如果需要修改这些分类, 你可以从库文件目录复制一份文件到你的工作目录, 然后编辑这个副本.

gmx select

目前, 一些分析工具支持一个扩展的概念 (动态)选择, 它与传统的索引组主要存在三个区别:

通过文本指定选择, 而不是从文件中读入固定的原子编号. 选择使用的语法类似于VMD. 文本可以交互地输入, 也可以过命令行提供, 或是使用文件.

选择并不局限于原子, 也可以是分析的对象, 如一组原子的质心位置. 一些工具可能并不支持那些非单个原子的选择, 例如, 当它们需要的信息, 如原子名称或类型之类, 只对单个原子可用时.

选择可以是动态的, 即, 对不同的轨迹帧选择不同的原子. 这样就可以只对体系满足一定的几何条件的那部分进行分析.

作为动态选择的一个简单例子, resname ABC and within 2 of resname DEF将选择残基名称为ABC, 且距残基DEF中任何一个原子都不超过2 nm的所有原子.

对支持动态选择的工具也可以使用传统的索引文件, 像以前的工具一样: 可以为工具指定一个.ndx文件, 通过组编号或组名称直接从索引文件中选择一个组. 索引组也可以作为更复杂选择的一部分使用.

作为开始, 你可以使用单个结构运行gmx select, 利用交互的提示试验不同的选择. 这个工具提供了一些输出选项, -on和-ofpdb可以分别将选择的原子写入索引文件和.pdb文件中. 这些选项不允许测试质心位置的选择, 但可以测试其他选择, 并检查结果.

如果想了解选择的详细语法和各个关键词, 对任何支持选择的工具你可以在交互提示下键入help, 或使用gmx help selections命令. 帮助被划分为不同的子主题, 每个子主题可以通过像help syntax / gmx help selections syntax 这样的命令查看. 一些单独的选择关键词还有其自己的扩展帮助, 可通过像help keywords within这样的命令进行查看.

目前的交互提示没有提供很多的编辑功能, 如果需要它们, 你可以在rlwrap下运行程序.

对那些不支持选择语法的工具, 你可以使用gmx select -on生成静态的索引组并传递给工具. 但是, 这种作法只能使用那些完全支持动态选择的工具所提供的灵活性中的一小部分(仅仅前面三项中的第一项).

为了利用这些选择的优点, 你也可以编写自己的分析工具, 具的例子请参看安装目录下share/gromacs/template中的template.cpp文件.

gmx view

在进行分析之前, 先查看一下你的轨迹通常很有帮助. GROMACS提供了一个简单的轨迹查看器gmx view. 这个查看器的优点是它不需要OpenGL的支持, 有些机器–如一些超级计算机–通常不支持OpenGL. 你也可以利用查看器生成EPS格式的轨迹(图 8.1). 如果你需要更快更花哨的查看器, 有一些程序可以读取GROMACS格式的轨迹, 请查看我们主页上更新的链接.

gmx energy, gmx traj

如果需要分析部分或全部的 能量 和其他性质, 如 总压力, 压力张量, 密度, 盒子体积 和 盒子尺寸, 请使用 gmx energy程序. 可以从一系列能量列表中进行选择, 如势能, 动能或总能, 也可以从对能量单独贡献的部分中进行选择, 如Lennard-Jones或二面角能量.

质心速度定义为

\[\bi V_{com} = {1\over M} \sum_{i=1}^N m_i \bi V_i \tag{8.2}\]

其中 \(M=\sum_{i=1}^N m_i\) 为体系的总质量. 可以利用程序gmx traj -com -ov对质心进行监测. 然而我们建议模拟时在每一步都移除质心的速度(参看第三章)!

gmx rdf

A类型粒子与B类型粒子之间的 径向分布函数 (RDF)或称对关联函数 \(g_{AB}(r)\) 的定义为

\[\alg g_{AB}(r) &= { \over {}_{local} } \\ &= { 1\over {}_{local} } {1\over N_A} \sum_{i \in A}^{N_A} \sum_{j \in B}^{N_B} {\delta (r_{ij}-r) \over 4\pi r^2} \tag{8.3} \ealg\]

其中, \(\) 为距离A粒子 \(r\) 处B粒子的密度, \(_{local}\) 为所有以A离子为中心, 半径 \(r_{max}\) 的壳层内B离子的平均密度(图 8.2C)

通常 \(r_{max}\) 的值取为盒子长度的一半. 计算时也会对时间进行平均. 具体计算时, 分析程序gmx rdf会将体系划分为球形切片(从 \(r\) 到 \(r+dr\), 图 8.2A), 并生成一个直方图而不是 \(\d\) 函数. SPC水模型氧-氧RDF[81]的一个例子如图 8.3.

也可以使用gmx rdf计算角度相关的RDF \(g_{AB}(r,\q)\), 其中角度 \(\q\) 根据特定的实验室坐标轴 \(\bi e\) 进行定义, 参看图8.2B.

\[\alg g_{AB}(r,\q) &= {1 \over {}_{local,\q} } {1 \over N_A} \sum_{i \in A}^{N_A} \sum_{j \in B}^{N_B} {\delta(r_{ij}-r)\delta(\theta_{ij}-\theta) \over 2\pi r^2 \sin(\theta)} \tag{8.4} \\ \cos(\theta_{ij}) &= {\bi r_{ij} \cdot \bi e \over \|r_{ij}\| \|e\|} \tag{8.5} \ealg\]

\(g_{AB}(r,\q)\) 在分析各向异性体系时很有用. 注意 在这种情况下, 归一化因子 \(_{local,\q}\) 为一直到 \(r_{max}\) 的所有从 \(\q\) 到 $\q+d\q$的角度切片的平均密度, 因此与角度相关. 参看图 8.2D

文献[106]对相关函数的理论作了很好的介绍. 我们在这里只介绍GROMACS代码中对相关函数实现的各种细节. 性质 \(f(t)\) 的自相关函数(ACF) \(C_f(t)\) 的定义为

\[C_f(t) = {< f(\xi)f(\xi+t) >}_\xi \tag{8.6}\]

其中, 等号右边的记号表示对 \(\x\) 进行平均, 即对不同时间起点进行平均. 对两种性质 \(f(t)\) 和 \(g(t)\) 也可以计算交叉相关函数

\[C_{fg}(t) = {< f(\xi)g(\xi+t) >}_\xi \tag{8.7}\]

然而, 在GROMACS中不存在计算相关函数的标准方法( 注意: 你可以使用xmgr程序计算交叉相关). 相关函数对时间的积分为相关时间 \(\t_f\):

\[\tau_f = \int_0^\infty C_f(t) \rmd t \tag{8.8}\]

具体计算时, 相关函数是根据时间间隔为 \(\D t\) 的离散数据的计算的, 因此, MD模拟的ACF为

\[C_f(j \D t) = {1 \over N-j} \sum_{i=0}^{N-1-j} f(i\D t) f((i+j)\D t) \tag{8.9}\]

其中, \(N\) 为计算时所有可用的时间帧. 所得的ACF明显地只适用于具有相同时间间隔 \(\D t\) 的时间点. 对许多应用, 由于需要知道ACF短时间的行为(例如, 开始10 ps), 这常常意味着我们保存数据的时间间隔必须比感兴趣的时间尺度短得多. 方程8.9隐含的另一点是, 理论上我们不可能以相同精度计算所有时间点上的ACF, 因为对 \(C_f(\D t)\) 我们有 \(N-1\) 个点, 而对 \(C_f((N-1)\D t)\) 只有1个点. 然而, 如果我们只计算长度为 \(M \D t\) 的ACF, 其中 \(M\le N/2\), 我们可以同样的统计精度计算所有的点

\[C_f(j\D t) = {1 \over M} \sum_{i=0}^{N-1-M} f(i\D t) f((i+j)\D t) \tag{8.10}\]

当然, 这里的 \(j \le M\). \(M\) 有时被视为相关函数的时间延迟. 当决定这样做时, 对非常短的时间间隔(\(j \ll M\)), 我们故意不使用所有的点, 但这使得我们更容易解释结果. 当从模拟轨迹计算ACF的时候, 不能忽略的另一点是, 通常时间起点 \(\x\) (方程8.6)并不是统计无关的, 这可能使得结果存在偏离. 这可以利用块平均方法进行测试.

在这种方法中, 计算时只包含间距最少为时间延迟的时间起点, 例如, 使用 \(k\) 乘上间距为 \(M\D t\) 的起点(其中 \(kM \le N\)):

\[C_f(j\D t) = {1 \over k}\sum_{i=0}^{k-1} f(i M\D t) f((iM+j) \D t) \tag{8.11}\]

然而, 利用这种方式, 需要非常长的模拟才能得到高精确度的结果, 因为对ACF有贡献的点的个数少得多.

根据方8.9计算ACF的代价正比于 \(N^2\), 计算量很可观. 然而, 通过使用快速Fourier变换计算卷积可以改进计算ACF的代价[106].

ACF存在一些重要的变形, 例如, 矢量 \(\bi p\) 的ACF

\[C_{\bi p}(t) = \int_0^\infty P_n(\cos \angle( \bi p(\xi), \bi p(\xi+t) ) \rmd \xi \tag{8.12}\]

其中 \(P_n(x)\) 为 \(n\) 阶Legendre多项式[1]. 这样的相关时间实际上可以NMR或其他弛豫实验中获得. GROMACS可以使用一阶和二阶Legendre多项式计算相关(方程8.12). 这种方法也可以用于旋转自相关(gmx rotacf)和偶极自相关(gmx dipoles).

为研究扭转角的动力学, 我们可以定义二面角的自相关函数[153]

\[C(t) = {}_{\tau} \tag{8.13}\]

注意 这不是两个函数的乘积, 像通常的相关函数那样, 但可以将它重写成两个乘积加和的形式:

\[C_(t) = {) M_{jj}^{1/2}(x_j-< x_j >) \right> \tag{8.22}\]

其中 \(M\) 为包含原子质量的对角矩阵(质量加权分析)或单位矩阵(非质量加权分析). \(C\) 是对称的 \(3N \times 3N\) 矩阵, 可以利用正交变换矩阵 \(R\) 对其进行对角化

\[R^T C R = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_{3N}),\qquad \lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_{3N} \tag{8.23}\]

\(R\) 的列就是本征矢量, 也称为主模或基本模. \(R\) 定义了一个到新坐标系的变换. 将轨迹投影到主模上便得到主成分 \(p_i(t)\)

\[\bi p(t) = R^T M^{1/2}(\bi x(t)-) \tag{8.24}\]

本征值 \(\l_i\) 为主成分 \(i\) 的均方涨落, 头几个主模通常描述了体系聚集的整体运动. 可以利用一个(或多个)主模对轨迹进行过滤. 对一个主模 \(i\), 其方式如下

\[\bi x^f(t)=< \bi x > + M^{-1/2} R_{*i} p_i(t) \tag{8.25}\]

当对大分子进行分析时, 通常想要移除整体的转动和平动以便仅仅观察内部运动. 这可以通过最小二乘拟合到参考结构实现. 需要注意的是, 参考结构对系综需要具有代表性, 因为参考结构的选择会影响协方差矩阵.

始终都应该检查是否已经很好地定义了主模. 如果第一个主成分类似于半个余弦, 第二个主成分类似于整个余弦, 你可能正在过滤噪声(参见下文). 检查头几个主模相关的一个好方式是计算模拟前半部分和后半部分采样的重叠. 注意 这只适用于两部分使用相同的参考结构的情况.

文献[159]定义了一个很好的重叠度量. 协方差矩阵的元素正比于位移的平方, 因此需要取矩阵的平方根来检查取样的范围. 可以从本征值 \(\l_i\) 和本征矢量计算矩阵的平方根, 本征矢量对应于矩阵 \(R\) 的列. 对一个对称, 对角占优, 大小为 \(3N \times 3N\) 的矩阵 \(A\), 其平方根为

\[A^{1/2} = R \text{diag}(\lambda_1^{1/2}, \lambda_2^{1/2}, \dots, \lambda_{3N}^{1/2}) R^T \tag{8.26}\]

容易证明这个矩阵及其自身的乘积为矩阵 \(A\). 现在我们可以定义协方差矩阵 \(A\) 和 \(B\) 之间的差异 \(d\)

\[\alg d(A,B) &= \sqrt{\text{tr} \left( (A^{1/2}-B^{1/2})^2 \right)} \tag{8.27} \\ &= \sqrt{\text{tr} \left( A+B-2A^{1/2}B^{1/2} \right)} \tag{8.28} \\ &= \left( \sum_{i=1}^N (\lambda_i^A + \lambda_i^B) - 2 \sum_{i=1}^N \sum_{j=1}^N \sqrt{\lambda_i^A \lambda_j^B} (R_i^A \cdot R_j^B)^2 \right)^{1/2} \tag{8.29} \ealg\]

其中, \(\text{tr}\) 为矩阵的迹. 重叠 \(s\) 的定义为

\[s(A,B) = 1 - {d(A,B) \over \sqrt{\text{tr} A + \text{tr} B} } \tag{8.30}\]

当且仅当矩阵 \(A\) 和 \(B\) 全等时, 重叠为1. 当抽样的子空间完全正交时, 重叠为0.

一个常用的度量是协方差矩阵头几个本征矢量子空间的重叠. \(m\) 个正交矢量 \(\bi w_1, \cdots, \bi w_m\) 张成的子空间和 \(n\) 个正交矢量 \(\bi v_1, \cdots, \bi v_n\) 张成的参考子空间之间的重叠可定量为

\[\text{overlap}(\bi v, \bi w) = {1\over n}\sum_{i=1}^n \sum_{j=1}^m (\bi v_i \cdot \bi w_j)^2 \tag{8.31}\]

重叠会随着 \(m\) 的增加而增大, 当集合 \(\bi v\) 为集合 \(\bi w\) 的子空间时, 重叠为1. 这个方法的缺点在于没有考虑到本征值. 所有本征矢量的权重都相同, 当存在简并的子空间时(本质值相等), 计算的重叠太小.

另一个有用的检查是余弦含量. 已经证明随机扩散的主成分是余弦, 且其周期数等于主成分索引数的一半[160,159]. 本征值正比于索引数的–2次方. 余弦含量定义为

\[{2 \over T}\left( \int_0^T \cos\left({i\pi t \over T}\right) p_i(t) \rmd t \right)^2 \left( \int_0^T p_i^2(t) \rmd t \right)^{-1} \tag{8.32}\]

当头几个主成分的余弦含量接近1时, 最大的涨落与势能没有联系, 而是对应于随机扩散.

gmx covar可以构建并对角化协方差矩阵. 可以利用gmx anaeig对主成分和重叠(还有许多其他量)进行绘图和分析. 余弦含量可以利用gmx analyze计算.

gmx angle, gmx covar, gmx anaeig

使用GROMACS可以在二面角空间[161]中进行主成分分析. 首先使用gmx mk_angndx或其他方法, 在索引文件中定义要计算的二面角, 然后使用gmx angle程序及其-or选项产生一个新的.trr文件, 其中包含了每个二面角分别在两个坐标系下的正弦和余弦值. 即, 在.trr文件中, 有相应于 \(\cos(\f_1)\), \(\sin(\f_1)\), \(\cos(\f_2)\), \(\sin(\f_2)\), \(\cdots\), \(\cos(\f_n)\), \(\sin(\f_n)\) 的一系列数字, 必要时, 会对数字序列补零. 然后, 你可以使用这个.trr文件作为gmx covar程序的输入, 像通常一样进行主成分分析. 为此, 你需要生成一个参考文件(.tpr,.gro,.pdb等), 其中包含与新的.trr文件相同数目的“原子”, 即, 对 \(n\) 个二面角, 你需要 \(2n/3\) 个原子(如果不是整数就舍去小数部分). 由于哑参考文件中的坐标与.trr文件中的信息毫无对应关系, 你需要使用gmx covar程序的-nofit选项. 对结果的分析可使用gmx anaeig.

gmx hbond

gmx hbond程序用于分析所有可能的施体D和受体A之间的 氢键 (HB). 分析时使用几何准则决定氢键的存在与否, 参看图 8.8:

\[\begin{split} r &\le r_{HB} = 0.35 \ \text{nm} \\ \alpha &\le \alpha_{HB} = 30^{\circ} \end{split} \tag{8.33}\]

参考值 \(r_{HB}=0.35 \text{nm}\) 对应于SPC水模型RDF的第一极小位置(参看图 8.3).

gmx hbond程序以下面的方式分析两组原子(它们必须相同或没有重叠)或指定的施体-氢-受体之间所有可能存在的氢键:

所有氢键中施体-受体距离的分布

所有氢键中氢-施体-受体之间夹角的分布

每个时间帧中氢键的总数目

残基之间实时的氢键数目, 划分到组 \(n-n+i\), 其中 \(n\) 和 \(n+i\) 代表残基编号, \(i\) 从0到6. \(i=6\) 的组也包括 \(i>6\) 的所有氢键. 这些组包括 \(n-n+3, n-n+4\) 和 \(n-n+5\) 氢键, 提供了 \(\a\) 螺旋, \(\b\) 转角或 \(\b\) 股形成的度量.

对所有氢键生存函数(0或1)的自相关函数进行平均, 可以计算氢键的寿命

\[C(\tau) = < s_i(t) s_i(t+\tau) > \tag{8.34}\]

对 \(t\) 时刻的氢键 \(i\), \(s_i(t)=\{0,1\}\). \(C(\t)\) 的积分给出了氢键平均寿命 \(\t_{HB}\) 的粗略估计

\[\tau_{HB} = \int_0^\infty C(\tau) \rmd \tau \tag{8.35}\]

程序会输出积分和完整的自相关函数 \(C(\t)\), 以便能够进行更复杂的分析(例如使用多指数拟合)得到 \(\t\) 的更佳估计值. 参考文献[162]给出了一个更完整的分析, 一个更花哨的选项是氢键动力学的Luzar和Chandler分析[163, 164].

生成大小为 #HB×#frames 的氢键存在映射图. 顺序与索引文件中的完全相同(见下文), 但反过来, 这意味着索引文件中的最后一个三联原子编号相应于存在映射图的第一行.

输出包含分析组的索引组, 这些组中所有的施体-氢原子对和受体原子, 分析组之间氢键中的施体-氢-受体三联对, 所有插入的溶剂分子.

gmx do_dssp, gmx rama, gmx wheel

为分析蛋白质结构的变化, 你可以计算回旋半径或随时间变化的最小残基距离(参看8.8节), 或计算RMSD(参看8.9节).

在运行中你也可以观察 二级结构元素 的变化, 为此, 你可以使用gmx do_dssp程序, 它是商业程序DSSP[165]的接口界面. 更多的信息请参看DSSP的手册. 图 8.10 给出了一个典型的gmx do_dssp输出图.

另一个重要的蛋白质分析是所谓的 Ramachandran图, 它是结构在蛋白质骨架的两个二面角 \(\f\) 和 \(\y\) 形成的平面上的投影. 参看图8.11.

你可以使用gmx rama获得Ramachandran图. 图 8.12是一个典型的输出.

当研究 \(\a\) 螺旋时, 多肽的 旋轮投影 很有用, 它可用于查看一个多肽是否是两性的. 你可以使用gmx wheel程序获得这样的投影图. 图 8.13给出了两个例子.

gmx order, gmx density, gmx potential, gmx traj

当模拟具有长碳尾的分子时, 可以计算它们的平均取向. 有许多不同的序参数, 它们大多数是相关的. gmx order程序可利用下面的方程计算序参数

\[S_z = {3 \over 2} -{1 \over 2} \tag{8.36}\]

其中 \(\q_z\) 为模拟盒子的 \(z\) 轴与待研究的分子轴之间的夹角. 后者被定义为从 \(\text C_{n-1}\) 到 \(\text C_{n+1}\) 的矢量. 可以同样的方式定义参数 \(S_x\) 和 \(S_y\). 尖括号隐含着要对时间和分子进行平均. 序参数可以在1(沿界面法向完全有序)和–1/2(完全垂直于界面法向)之间变化, 值为零时代表取向具有各向同性.

程序可以完成两件事. 它可以分别计算每个 \(\ce{CH2}\) 片段的序参数, 对所有三个轴的序参数, 或将盒子划分为切片并计算一个切片内每个片段序参数的平均值.

第一个方法给出了分子从头到尾如何排列的说明, 第二个方法给出了排列和盒子长度的关系.

从一条轨迹中可以计算出穿过界面的静电势(\(\y\)), 通过电荷密度(\(\r(z)\))的双重积分

\[\psi(z)-\psi(-\infty) = - \int_{-\infty}^z \rmd z' \int_{-\infty}^{z'} \r(z'') \rmd z''/\e_0 \tag{8.37}\]

其中, 位置 \(z =-\infty\) 在体相中足够远以致场强为零. 使用这种方法, 可以将总的势能“划分”为层和水分子的独立贡献. gmx potential程序将盒子划分为切片, 并对每个切片中原子的所有电荷进行加和. 然后程序会积分电荷密度得到电场, 对电场进行积分给出势能. 电荷密度, 电场和势能会写入xvgr输入文件中.

gmx traj是一个非常简单的分析程序, 它所做的是打印选中原子的坐标, 速度, 或力. 它也可以计算一个或多个分子的质心, 并将质心坐标打印到三个文件中. 这个分析程序本身可能并不是很有用, 但获得了选择的分子或原子的坐标后很利于进行进一步的分析, 不仅仅对于界面体系是这样.

gmx density程序计算组的质量密度, 并给出沿盒子轴向的密度分布. 对于观察组的分布或穿过界面的原子, 这很有用.

\(P_0(x)=1, P_1(x)=x, P_2(x)=(3x^2-1)/2\)  ↩